TPC fitting
20220829
PR data
Import data
data <- read.csv("output/pr_curve_extracted_rates.csv")
str(data)## 'data.frame': 1114 obs. of 9 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Species : chr "P. grandis" "P. grandis" "P. grandis" "P. grandis" ...
## $ colony_id.x : chr "POG-150" "POG-150" "POG-150" "POG-150" ...
## $ avg_temp_interval: num 17.2 17.1 21.1 21.1 24.1 ...
## $ Temperature : int 17 17 21 21 24 24 28 28 30 30 ...
## $ Light : int 0 560 0 560 0 560 0 560 0 560 ...
## $ Run : int 4 4 4 4 4 4 4 4 4 4 ...
## $ micromol.cm2.s : num -0.0583 0.0334 -0.0762 0.0421 -0.0981 ...
## $ micromol.cm2.h : num -210 120 -274 152 -353 ...data$Temperature <- as.numeric(data$avg_temp_interval)
str(data)## 'data.frame': 1114 obs. of 9 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Species : chr "P. grandis" "P. grandis" "P. grandis" "P. grandis" ...
## $ colony_id.x : chr "POG-150" "POG-150" "POG-150" "POG-150" ...
## $ avg_temp_interval: num 17.2 17.1 21.1 21.1 24.1 ...
## $ Temperature : num 17.2 17.1 21.1 21.1 24.1 ...
## $ Light : int 0 560 0 560 0 560 0 560 0 560 ...
## $ Run : int 4 4 4 4 4 4 4 4 4 4 ...
## $ micromol.cm2.s : num -0.0583 0.0334 -0.0762 0.0421 -0.0981 ...
## $ micromol.cm2.h : num -210 120 -274 152 -353 ...data <- data %>%
arrange(colony_id.x) %>%
group_by(colony_id.x) %>%
mutate(curve_id = group_indices())
data %>%
ggplot(aes(x=Temperature, y=micromol.cm2.s, color=Species, group=colony_id.x))+
geom_point()+
geom_line()+
ylab("µmol O2 cm^2 s-1")+
facet_wrap(~ Species*Light, ncol = 5)+
theme(legend.position = "none")
d <- data %>%
filter(Light==0) %>%
filter(!Species == "Recruit")
d$micromol.cm2.s <- -d$micromol.cm2.s
d$micromol.cm2.s <- replace(d$micromol.cm2.s, d$micromol.cm2.s<0,0)
d <- d[,-c(1,4,6,7,9)]
str(d)## gropd_df [278 × 5] (S3: grouped_df/tbl_df/tbl/data.frame)
## $ Species : chr [1:278] "P. grandis" "P. grandis" "P. grandis" "P. grandis" ...
## $ colony_id.x : chr [1:278] "POG-150" "POG-150" "POG-150" "POG-150" ...
## $ Temperature : num [1:278] 17.2 21.1 24.1 28 29.7 ...
## $ micromol.cm2.s: num [1:278] 0.0583 0.0762 0.0981 0.1325 0.1434 ...
## $ curve_id : int [1:278] 1 1 1 1 1 1 1 1 1 2 ...
## - attr(*, "groups")= tibble [30 × 2] (S3: tbl_df/tbl/data.frame)
## ..$ colony_id.x: chr [1:30] "POG-150" "POG-151" "POG-156" "POG-159" ...
## ..$ .rows : list<int> [1:30]
## .. ..$ : int [1:9] 1 2 3 4 5 6 7 8 9
## .. ..$ : int [1:9] 10 11 12 13 14 15 16 17 18
## .. ..$ : int [1:8] 19 20 21 22 23 24 25 26
## .. ..$ : int [1:8] 27 28 29 30 31 32 33 34
## .. ..$ : int [1:9] 35 36 37 38 39 40 41 42 43
## .. ..$ : int [1:8] 44 45 46 47 48 49 50 51
## .. ..$ : int [1:9] 52 53 54 55 56 57 58 59 60
## .. ..$ : int [1:9] 61 62 63 64 65 66 67 68 69
## .. ..$ : int [1:9] 70 71 72 73 74 75 76 77 78
## .. ..$ : int [1:8] 79 80 81 82 83 84 85 86
## .. ..$ : int [1:8] 87 88 89 90 91 92 93 94
## .. ..$ : int [1:8] 95 96 97 98 99 100 101 102
## .. ..$ : int [1:8] 103 104 105 106 107 108 109 110
## .. ..$ : int [1:9] 111 112 113 114 115 116 117 118 119
## .. ..$ : int [1:8] 120 121 122 123 124 125 126 127
## .. ..$ : int [1:9] 128 129 130 131 132 133 134 135 136
## .. ..$ : int [1:9] 137 138 139 140 141 142 143 144 145
## .. ..$ : int [1:8] 146 147 148 149 150 151 152 153
## .. ..$ : int [1:8] 154 155 156 157 158 159 160 161
## .. ..$ : int [1:9] 162 163 164 165 166 167 168 169 170
## .. ..$ : int [1:9] 171 172 173 174 175 176 177 178 179
## .. ..$ : int [1:9] 180 181 182 183 184 185 186 187 188
## .. ..$ : int [1:7] 189 190 191 192 193 194 195
## .. ..$ : int [1:7] 196 197 198 199 200 201 202
## .. ..$ : int [1:32] 203 204 205 206 207 208 209 210 211 212 ...
## .. ..$ : int [1:9] 235 236 237 238 239 240 241 242 243
## .. ..$ : int [1:9] 244 245 246 247 248 249 250 251 252
## .. ..$ : int [1:9] 253 254 255 256 257 258 259 260 261
## .. ..$ : int [1:8] 262 263 264 265 266 267 268 269
## .. ..$ : int [1:9] 270 271 272 273 274 275 276 277 278
## .. ..@ ptype: int(0)
## ..- attr(*, ".drop")= logi TRUEcolnames(d) <- c("Species", "colony_id", "temp","rate","curve_id")
d$curve_id <- as.numeric(d$curve_id)
str(d)## gropd_df [278 × 5] (S3: grouped_df/tbl_df/tbl/data.frame)
## $ Species : chr [1:278] "P. grandis" "P. grandis" "P. grandis" "P. grandis" ...
## $ colony_id: chr [1:278] "POG-150" "POG-150" "POG-150" "POG-150" ...
## $ temp : num [1:278] 17.2 21.1 24.1 28 29.7 ...
## $ rate : num [1:278] 0.0583 0.0762 0.0981 0.1325 0.1434 ...
## $ curve_id : num [1:278] 1 1 1 1 1 1 1 1 1 2 ...
## - attr(*, "groups")= tibble [30 × 2] (S3: tbl_df/tbl/data.frame)
## ..$ colony_id: chr [1:30] "POG-150" "POG-151" "POG-156" "POG-159" ...
## ..$ .rows : list<int> [1:30]
## .. ..$ : int [1:9] 1 2 3 4 5 6 7 8 9
## .. ..$ : int [1:9] 10 11 12 13 14 15 16 17 18
## .. ..$ : int [1:8] 19 20 21 22 23 24 25 26
## .. ..$ : int [1:8] 27 28 29 30 31 32 33 34
## .. ..$ : int [1:9] 35 36 37 38 39 40 41 42 43
## .. ..$ : int [1:8] 44 45 46 47 48 49 50 51
## .. ..$ : int [1:9] 52 53 54 55 56 57 58 59 60
## .. ..$ : int [1:9] 61 62 63 64 65 66 67 68 69
## .. ..$ : int [1:9] 70 71 72 73 74 75 76 77 78
## .. ..$ : int [1:8] 79 80 81 82 83 84 85 86
## .. ..$ : int [1:8] 87 88 89 90 91 92 93 94
## .. ..$ : int [1:8] 95 96 97 98 99 100 101 102
## .. ..$ : int [1:8] 103 104 105 106 107 108 109 110
## .. ..$ : int [1:9] 111 112 113 114 115 116 117 118 119
## .. ..$ : int [1:8] 120 121 122 123 124 125 126 127
## .. ..$ : int [1:9] 128 129 130 131 132 133 134 135 136
## .. ..$ : int [1:9] 137 138 139 140 141 142 143 144 145
## .. ..$ : int [1:8] 146 147 148 149 150 151 152 153
## .. ..$ : int [1:8] 154 155 156 157 158 159 160 161
## .. ..$ : int [1:9] 162 163 164 165 166 167 168 169 170
## .. ..$ : int [1:9] 171 172 173 174 175 176 177 178 179
## .. ..$ : int [1:9] 180 181 182 183 184 185 186 187 188
## .. ..$ : int [1:7] 189 190 191 192 193 194 195
## .. ..$ : int [1:7] 196 197 198 199 200 201 202
## .. ..$ : int [1:32] 203 204 205 206 207 208 209 210 211 212 ...
## .. ..$ : int [1:9] 235 236 237 238 239 240 241 242 243
## .. ..$ : int [1:9] 244 245 246 247 248 249 250 251 252
## .. ..$ : int [1:9] 253 254 255 256 257 258 259 260 261
## .. ..$ : int [1:8] 262 263 264 265 266 267 268 269
## .. ..$ : int [1:9] 270 271 272 273 274 275 276 277 278
## .. ..@ ptype: int(0)
## ..- attr(*, ".drop")= logi TRUETPC fitting Padifeld et al rTPC and nls.multstart: A new
pipeline to fit thermal performance curves in r
https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13585
Sharpe Schoolfield 1981 model Schoolfield, R. M., Sharpe, P. J. H., & Magnuson, C. E. (1981). Non-linear regression of biological temperature-dependent rate models based on absolute reaction-rate theory. Journal of theoretical biology, 88(4), 719-731. https://doi.org/10.1016/0022-5193(81)90246-0
#edit nls_multstart to allow for a progress bar
# edit nls_multstart to allow for a progress bar
nls_multstart_progress <- function(formula, data = parent.frame(), iter, start_lower,
start_upper, supp_errors = c("Y", "N"), convergence_count = 100,
control, modelweights, ...){
if(!is.null(pb)){
pb$tick()
}
nls_multstart(formula = formula, data = data, iter = iter, start_lower = start_lower,
start_upper = start_upper, supp_errors = supp_errors, convergence_count = convergence_count,
control = control, modelweights = modelweights, ...)
}start progress bar and estimate time it will take
# start progress bar and estimate time it will take
number_of_models <- 1
number_of_curves <- length(unique(d$curve_id))
# setup progress bar
pb <- progress::progress_bar$new(total = number_of_curves*number_of_models,
clear = FALSE,
format ="[:bar] :percent :elapsedfull")fit chosen model formulation in rTPC
# fit chosen model formulation in rTPC
d_fits <- nest(d, data = c(temp, rate)) %>%
mutate(sharpeschoolhigh = map(data, ~nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = .x,
iter = c(3,3,3,3),
start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'sharpeschoolhigh_1981') - 1,
start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'sharpeschoolhigh_1981') + 1,
lower = get_lower_lims(.x$temp, .x$rate, model_name = 'sharpeschoolhigh_1981'),
upper = get_upper_lims(.x$temp, .x$rate, model_name = 'sharpeschoolhigh_1981'),
supp_errors = 'Y',
convergence_count = FALSE)))d_preds <- mutate(d_fits, new_data = map(data, ~tibble(temp = seq(min(.x$temp), max(.x$temp), length.out = 100)))) %>%
# get rid of original data column
select(., -data) %>%
# stack models into a single column, with an id column for model_name
pivot_longer(., names_to = 'model_name', values_to = 'fit', c(sharpeschoolhigh)) %>%
# create new list column containing the predictions
# this uses both fit and new_data list columns
mutate(preds = map2(fit, new_data, ~augment(.x, newdata = .y))) %>%
# select only the columns we want to keep
select(curve_id, Species,model_name, preds) %>%
# unlist the preds list column
unnest(preds)
glimpse(d_preds)## Rows: 3,000
## Columns: 6
## Groups: colony_id [30]
## $ colony_id <chr> "POG-150", "POG-150", "POG-150", "POG-150", "POG-150", "POG…
## $ curve_id <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ Species <chr> "P. grandis", "P. grandis", "P. grandis", "P. grandis", "P.…
## $ model_name <chr> "sharpeschoolhigh", "sharpeschoolhigh", "sharpeschoolhigh",…
## $ temp <dbl> 17.21817, 17.44516, 17.67214, 17.89913, 18.12611, 18.35309,…
## $ .fitted <dbl> 0.05791013, 0.05903478, 0.06017729, 0.06133772, 0.06251611,…ggplot(d_preds) +
geom_line(aes(temp, .fitted, col = Species, group=curve_id)) +
geom_point(aes(temp, rate, col = Species), d) +
#facet_wrap(~curve_id*Species, scales = 'free_y', ncol = 6) +
theme_bw() +
#theme(legend.position = 'none') +
scale_color_brewer(type = 'qual', palette = 2) +
labs(x = 'Temperature (ºC)',
y = 'Light Enhanced Dark Respiration rate',
title = 'sharpeschoolhigh')
ggplot(d_preds) +
geom_line(aes(temp, .fitted, col = Species, group=curve_id)) +
geom_point(aes(temp, rate, col = Species), d) +
facet_wrap(~Species, ncol = 4) +
theme_bw() +
theme(legend.position = 'none') +
scale_color_brewer(type = 'qual', palette = 2) +
labs(x = 'Temperature (ºC)',
y = 'Light Enhanced Dark Respiration rate',
title = 'sharpeschoolhigh')
#calculate TPC metrics by species
d_params <- pivot_longer(d_fits, names_to = 'model_name', values_to = 'fit', c(sharpeschoolhigh)) %>%
mutate(params = map(fit, calc_params)) %>%
select(curve_id, Species, model_name, params) %>%
unnest(params)
glimpse(d_params)## Rows: 30
## Columns: 15
## Groups: colony_id [30]
## $ colony_id <chr> "POG-150", "POG-151", "POG-156", "POG-159", "POG…
## $ curve_id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1…
## $ Species <chr> "P. grandis", "P. grandis", "P. grandis", "P. gr…
## $ model_name <chr> "sharpeschoolhigh", "sharpeschoolhigh", "sharpes…
## $ rmax <dbl> 0.1341227, 0.1692973, 0.1834050, 0.1645734, 0.14…
## $ topt <dbl> 31.46017, 32.99975, 33.07686, 32.19303, 33.44543…
## $ ctmin <dbl> -5.66082708, -5.93625000, -5.69314000, -9.320970…
## $ ctmax <dbl> 55.55617, 48.59375, 44.51986, 54.66203, 47.32543…
## $ e <dbl> 0.5622984, 0.5017270, 0.4942046, 0.3974826, 0.48…
## $ eh <dbl> 0.3230061, 0.9464816, 1.6929276, NA, 0.9786143, …
## $ q10 <dbl> 2.088531, 1.917583, 1.890915, 1.675354, 1.883673…
## $ thermal_safety_margin <dbl> 24.096, 15.594, 11.443, 22.469, 13.880, 10.085, …
## $ thermal_tolerance <dbl> 61.217, 54.530, 50.213, 63.983, 54.726, 48.264, …
## $ breadth <dbl> 11.432, 9.517, 8.169, 11.732, 9.195, 7.654, 8.30…
## $ skewness <dbl> 0.23929229, -0.44475468, -1.19872297, NA, -0.490…TPC.pars <- as.data.frame(d_params)
summary(aov(rmax ~ Species, TPC.pars))## Df Sum Sq Mean Sq F value Pr(>F)
## Species 3 0.001081 0.0003603 0.868 0.47
## Residuals 26 0.010789 0.0004150summary(aov(ctmax ~ Species, TPC.pars))## Df Sum Sq Mean Sq F value Pr(>F)
## Species 3 98.2 32.73 1.772 0.177
## Residuals 26 480.2 18.47summary(aov(e ~ Species, TPC.pars))## Df Sum Sq Mean Sq F value Pr(>F)
## Species 3 0.0021 0.00070 0.05 0.985
## Residuals 26 0.3612 0.01389summary(aov(eh ~ Species, TPC.pars))## Df Sum Sq Mean Sq F value Pr(>F)
## Species 3 0.308 0.1028 0.296 0.828
## Residuals 17 5.910 0.3476
## 9 observations deleted due to missingnesssummary(aov(breadth ~ Species, TPC.pars))## Df Sum Sq Mean Sq F value Pr(>F)
## Species 3 9.92 3.306 1.888 0.156
## Residuals 26 45.52 1.751summary(aov(topt ~ Species, TPC.pars))## Df Sum Sq Mean Sq F value Pr(>F)
## Species 3 3.11 1.037 0.623 0.606
## Residuals 26 43.26 1.664d_fits <- mutate(d_fits, bootstrap = list(rep(NA, n())))
# run for loop to bootstrap each refitted model
for(i in 1:nrow(d_fits)){
temp_data <- d_fits$data[[i]]
temp_fit <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = temp_data,
start = coef(d_fits$sharpeschoolhigh[[i]]),
lower = get_lower_lims(temp_data$temp, temp_data$rate, model_name = 'sharpeschoolhigh_1981')-1,
upper = get_upper_lims(temp_data$temp, temp_data$rate, model_name = 'sharpeschoolhigh_1981')+1)
boot <- Boot(temp_fit, method = 'residual')
d_fits$bootstrap[[i]] <- boot
rm(list = c('temp_fit', 'temp_data', 'boot'))
}# get the raw values of each bootstrap
d_fits <- mutate(d_fits, output_boot = map(bootstrap, function(x) x$t))
# calculate predictions with a gnarly written function
d_fits <- mutate(d_fits, preds = map2(output_boot, data, function(x, y){
temp <- as.data.frame(x) %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(y$temp), max(y$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = sharpeschoolhigh_1981(temp, r_tref,e,eh,th, tref=28))
return(temp)
}))
# select, unnest and calculate 95% CIs of predictions
boot_conf_preds <- select(d_fits, curve_id, preds) %>%
unnest(preds) %>%
group_by(curve_id, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975),
.groups = 'drop')
ggplot() +
geom_line(aes(temp, .fitted), d_preds, col = 'blue') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), boot_conf_preds, fill = 'blue', alpha = 0.3) +
geom_point(aes(temp, rate), d, size = 2) +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Rate') +
facet_wrap(~curve_id, ncol=4)
# GRANDIS
# grandis CURVE FIT
# P. grandis
d.grandis <- d %>%
filter(Species == "P. grandis")
#fit
grandis.fit <- nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d,
iter = c(3,3,3,3),
start_lower = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') - 1,
start_upper = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') + 1,
lower = get_lower_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
upper = get_upper_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
supp_errors = 'Y',
convergence_count = FALSE)
grandis.fit## Nonlinear regression model
## model: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 28)
## data: data
## r_tref e eh th
## 0.1413 0.5775 3.4506 37.0533
## residual sum-of-squares: 0.1234
##
## Number of iterations to convergence: 14
## Achieved convergence tolerance: 1.49e-08#generate the predicted data
grandis_new_data <- data.frame(temp = seq(min(d.grandis$temp), max(d.grandis$temp), 0.5))
grandis.preds <- augment(grandis.fit, newdata = grandis_new_data)
#calculate TPC parameters
grandis.TCP.res <- calc_params(grandis.fit) %>%
mutate_all(round, 2) # round
grandis.TCP.res## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.17 33.24 -5.65 47.29 0.51 1.11 1.93 14.04 52.94
## breadth skewness
## 1 9 -0.6### Bootstrapping curve fit
# refit model using nlsLM
grandis.fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d.grandis,
start = coef(grandis.fit),
lower = get_lower_lims(d.grandis$temp, d.grandis$rate, model_name = 'sharpeschoolhigh_1981')-1,
upper = get_upper_lims(d.grandis$temp, d.grandis$rate, model_name = 'sharpeschoolhigh_1981')+1,
weights = rep(1, times = nrow(d.grandis)))
# bootstrap using case resampling
grandis.boot1 <- Boot(grandis.fit_nlsLM, method = 'case')
# look at the data
head(grandis.boot1$t)## r_tref e eh th
## [1,] 0.1305408 0.4820667 3.920797 38.09578
## [2,] 0.1417562 0.5525771 3.148593 37.36422
## [3,] 0.1389140 0.5897576 3.334068 37.19541
## [4,] 0.1334456 0.4320001 4.887677 38.02246
## [5,] 0.1437546 0.6846489 3.129592 36.34727
## [6,] 0.1450414 0.5766447 2.596644 36.64574# create predictions of each bootstrapped model
grandis.boot1_preds <- grandis.boot1$t %>%
as.data.frame() %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(d.grandis$temp), max(d.grandis$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 28))
# calculate bootstrapped confidence intervals
grandis.boot1_conf_preds <- group_by(grandis.boot1_preds, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975)) %>%
ungroup()
# plot bootstrapped CIs
grandis.CI.plot <- ggplot() +
geom_line(aes(temp, .fitted), grandis.preds, col = 'darkgreen') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), grandis.boot1_conf_preds, fill = 'darkgreen', alpha = 0.3) +
geom_point(aes(temp, rate), d.grandis, size = 2, alpha = 0.5,col = 'darkgreen') +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Rate (µmol O2/cm2/s)')
grandis.CI.plot
P. meandrina
# meandrina CURVE FIT
# P. meandrina
d.meandrina <- d %>%
filter(Species == "P. meandrina")
#fit
meandrina.fit <- nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d,
iter = c(3,3,3,3),
start_lower = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') - 1,
start_upper = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') + 1,
lower = get_lower_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
upper = get_upper_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
supp_errors = 'Y',
convergence_count = FALSE)
meandrina.fit## Nonlinear regression model
## model: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 28)
## data: data
## r_tref e eh th
## 0.1413 0.5775 3.4506 37.0533
## residual sum-of-squares: 0.1234
##
## Number of iterations to convergence: 14
## Achieved convergence tolerance: 1.49e-08#generate the predicted data
meandrina_new_data <- data.frame(temp = seq(min(d.meandrina$temp), max(d.meandrina$temp), 0.5))
meandrina.preds <- augment(meandrina.fit, newdata = meandrina_new_data)
#calculate TPC parameters
meandrina.TCP.res <- calc_params(meandrina.fit) %>%
mutate_all(round, 2) # round
meandrina.TCP.res## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.17 33.24 -5.65 47.29 0.51 1.11 1.93 14.04 52.94
## breadth skewness
## 1 9 -0.6### Bootstrapping curve fit
# refit model using nlsLM
meandrina.fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d.meandrina,
start = coef(meandrina.fit),
lower = get_lower_lims(d.meandrina$temp, d.meandrina$rate, model_name = 'sharpeschoolhigh_1981')-1,
upper = get_upper_lims(d.meandrina$temp, d.meandrina$rate, model_name = 'sharpeschoolhigh_1981')+1,
weights = rep(1, times = nrow(d.meandrina)))
# bootstrap using case resampling
meandrina.boot1 <- Boot(meandrina.fit_nlsLM, method = 'case')
# look at the data
head(meandrina.boot1$t)## r_tref e eh th
## [1,] 0.1509585 0.6786125 3.527396 36.04762
## [2,] 0.1375984 0.5540254 3.779009 37.61593
## [3,] 0.1431038 0.6306717 3.714936 36.51225
## [4,] 0.1456467 0.5567787 3.469522 36.88690
## [5,] 0.1316263 0.3965938 5.982946 38.59573
## [6,] 0.1478011 0.5880170 3.168019 36.62537# create predictions of each bootstrapped model
meandrina.boot1_preds <- meandrina.boot1$t %>%
as.data.frame() %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(d.meandrina$temp), max(d.meandrina$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 28))
# calculate bootstrapped confidence intervals
meandrina.boot1_conf_preds <- group_by(meandrina.boot1_preds, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975)) %>%
ungroup()
# plot bootstrapped CIs
meandrina.CI.plot <- ggplot() +
geom_line(aes(temp, .fitted), meandrina.preds, col = 'orange') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), meandrina.boot1_conf_preds, fill = 'orange', alpha = 0.3) +
geom_point(aes(temp, rate), d.meandrina, size = 2, alpha = 0.5,col = 'orange') +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Rate (µmol O2/cm2/s)')
meandrina.CI.plot
P. tuahiniensis
# tuahiniensis CURVE FIT
# P. tuahiniensis
d.tuahiniensis <- d %>%
filter(Species == "P. tuahiniensis")
#fit
tuahiniensis.fit <- nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d,
iter = c(3,3,3,3),
start_lower = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') - 1,
start_upper = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') + 1,
lower = get_lower_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
upper = get_upper_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
supp_errors = 'Y',
convergence_count = FALSE)
tuahiniensis.fit## Nonlinear regression model
## model: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 28)
## data: data
## r_tref e eh th
## 0.1413 0.5775 3.4506 37.0533
## residual sum-of-squares: 0.1234
##
## Number of iterations to convergence: 14
## Achieved convergence tolerance: 1.49e-08#generate the predicted data
tuahiniensis_new_data <- data.frame(temp = seq(min(d.tuahiniensis$temp), max(d.tuahiniensis$temp), 0.5))
tuahiniensis.preds <- augment(tuahiniensis.fit, newdata = tuahiniensis_new_data)
#calculate TPC parameters
tuahiniensis.TCP.res <- calc_params(tuahiniensis.fit) %>%
mutate_all(round, 2) # round
tuahiniensis.TCP.res## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.17 33.24 -5.65 47.29 0.51 1.11 1.93 14.04 52.94
## breadth skewness
## 1 9 -0.6### Bootstrapping curve fit
# refit model using nlsLM
tuahiniensis.fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d.tuahiniensis,
start = coef(tuahiniensis.fit),
lower = get_lower_lims(d.tuahiniensis$temp, d.tuahiniensis$rate, model_name = 'sharpeschoolhigh_1981')-1,
upper = get_upper_lims(d.tuahiniensis$temp, d.tuahiniensis$rate, model_name = 'sharpeschoolhigh_1981')+1,
weights = rep(1, times = nrow(d.tuahiniensis)))
# bootstrap using case resampling
tuahiniensis.boot1 <- Boot(tuahiniensis.fit_nlsLM, method = 'case')
# look at the data
head(tuahiniensis.boot1$t)## r_tref e eh th
## [1,] 0.1301138 0.5511504 4.414331 37.26193
## [2,] 0.1420862 0.6200535 4.002197 36.88074
## [3,] 0.1225529 0.5349330 5.719269 38.16553
## [4,] 0.1357566 0.5963970 4.068060 37.37621
## [5,] 0.1295459 0.5457687 5.117105 37.60311
## [6,] 0.1365303 0.6017160 3.325320 36.68036# create predictions of each bootstrapped model
tuahiniensis.boot1_preds <- tuahiniensis.boot1$t %>%
as.data.frame() %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(d.tuahiniensis$temp), max(d.tuahiniensis$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 28))
# calculate bootstrapped confidence intervals
tuahiniensis.boot1_conf_preds <- group_by(tuahiniensis.boot1_preds, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975)) %>%
ungroup()
# plot bootstrapped CIs
tuahiniensis.CI.plot <- ggplot() +
geom_line(aes(temp, .fitted), tuahiniensis.preds, col = 'purple') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), tuahiniensis.boot1_conf_preds, fill = 'purple', alpha = 0.3) +
geom_point(aes(temp, rate), d.tuahiniensis, size = 2, alpha = 0.5,col = 'purple') +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Rate (µmol O2/cm2/s)')
tuahiniensis.CI.plot
P. verrucosa
# verrucosa CURVE FIT
# P. verrucosa
d.verrucosa <- d %>%
filter(Species == "P. verrucosa")
#fit
verrucosa.fit <- nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d,
iter = c(3,3,3,3),
start_lower = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') - 1,
start_upper = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') + 1,
lower = get_lower_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
upper = get_upper_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
supp_errors = 'Y',
convergence_count = FALSE)
verrucosa.fit## Nonlinear regression model
## model: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 28)
## data: data
## r_tref e eh th
## 0.1413 0.5775 3.4506 37.0533
## residual sum-of-squares: 0.1234
##
## Number of iterations to convergence: 14
## Achieved convergence tolerance: 1.49e-08#generate the predicted data
verrucosa_new_data <- data.frame(temp = seq(min(d.verrucosa$temp), max(d.verrucosa$temp), 0.5))
verrucosa.preds <- augment(verrucosa.fit, newdata = verrucosa_new_data)
#calculate TPC parameters
verrucosa.TCP.res <- calc_params(verrucosa.fit) %>%
mutate_all(round, 2) # round
verrucosa.TCP.res## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.17 33.24 -5.65 47.29 0.51 1.11 1.93 14.04 52.94
## breadth skewness
## 1 9 -0.6### Bootstrapping curve fit
# refit model using nlsLM
verrucosa.fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d.verrucosa,
start = coef(verrucosa.fit),
lower = get_lower_lims(d.verrucosa$temp, d.verrucosa$rate, model_name = 'sharpeschoolhigh_1981')-1,
upper = get_upper_lims(d.verrucosa$temp, d.verrucosa$rate, model_name = 'sharpeschoolhigh_1981')+1,
weights = rep(1, times = nrow(d.verrucosa)))
# bootstrap using case resampling
verrucosa.boot1 <- Boot(verrucosa.fit_nlsLM, method = 'case')
# look at the data
head(verrucosa.boot1$t)## r_tref e eh th
## [1,] 0.1525225 0.5886983 3.320467 36.52780
## [2,] 0.1511853 0.5745659 3.419156 36.81460
## [3,] 0.1487430 0.5655326 3.291840 37.15557
## [4,] 0.1507273 0.5694067 3.428851 37.05959
## [5,] 0.1489431 0.5968072 3.058033 36.80796
## [6,] 0.1438597 0.5573666 4.150847 37.09051# create predictions of each bootstrapped model
verrucosa.boot1_preds <- verrucosa.boot1$t %>%
as.data.frame() %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(d.verrucosa$temp), max(d.verrucosa$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 28))
# calculate bootstrapped confidence intervals
verrucosa.boot1_conf_preds <- group_by(verrucosa.boot1_preds, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975)) %>%
ungroup()
# plot bootstrapped CIs
verrucosa.CI.plot <- ggplot() +
geom_line(aes(temp, .fitted), verrucosa.preds, col = 'pink') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), verrucosa.boot1_conf_preds, fill = 'pink', alpha = 0.3) +
geom_point(aes(temp, rate), d.verrucosa, size = 2, alpha = 0.5,col = 'pink') +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Rate (µmol O2/cm2/s)')
verrucosa.CI.plot
#set plot colors
cols <- c("grandis"="darkgreen","meandrina"="orange", "tuahiniensis"="purple", "verrucosa"="pink")
# plot data and model fit
TPC.plot <- ggplot(data=d, aes(x=temp)) +
geom_point(aes(temp, rate, color="grandis"), d.grandis, size = 2, alpha = 0.5) +
geom_point(aes(temp, rate, color="meandrina"), d.meandrina, size = 2, alpha = 0.5) +
geom_point(aes(temp, rate, color="tuahiniensis"), d.tuahiniensis, size = 2, alpha = 0.5) +
geom_point(aes(temp, rate, color="verrucosa"), d.verrucosa, size = 2, alpha = 0.5) +
geom_line(aes(temp, .fitted), grandis.preds, col = 'darkgreen', size=2) +
geom_line(aes(temp, .fitted), meandrina.preds, col = 'orange', size=2) +
geom_line(aes(temp, .fitted), tuahiniensis.preds, col = "purple", size=2) +
geom_line(aes(temp, .fitted), verrucosa.preds, col = "pink", size=2) +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), grandis.boot1_conf_preds, fill = "darkgreen", alpha = 0.3) +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), meandrina.boot1_conf_preds, fill = "orange", alpha = 0.3) +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), tuahiniensis.boot1_conf_preds, fill = 'purple', alpha = 0.3) +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), verrucosa.boot1_conf_preds, fill = 'pink', alpha = 0.3) +
xlim(16,41)+
#ylim(0,0.45)+
scale_x_continuous(breaks=c(16,22,24,26,28,30,32,34,36,41))+
theme_bw(base_size = 12) +
scale_colour_manual(name="Species",values=cols)+
theme(legend.position = "none",
panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
labs(x = 'Temperature (ºC)',
y = expression("Rate"~µmol~O[2] ~cm^{-2}~s^{-1}))
TPC.plot 
Bootstrap CI for all TPC parameters
broom::tidy(grandis.fit_nlsLM)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 r_tref 0.139 0.00609 22.8 1.16e-32
## 2 e 0.602 0.0742 8.12 1.82e-11
## 3 eh 3.18 0.464 6.86 3.07e- 9
## 4 th 37.0 0.726 51.0 3.91e-54broom::tidy(meandrina.fit_nlsLM)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 r_tref 0.141 0.00573 24.6 4.91e-31
## 2 e 0.566 0.0753 7.51 6.19e-10
## 3 eh 3.77 0.632 5.97 1.93e- 7
## 4 th 37.1 0.621 59.7 5.39e-51broom::tidy(tuahiniensis.fit_nlsLM)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 r_tref 0.132 0.00535 24.7 2.10e-34
## 2 e 0.562 0.0778 7.23 7.41e-10
## 3 eh 3.91 0.698 5.60 4.88e- 7
## 4 th 37.5 0.613 61.1 1.71e-58broom::tidy(verrucosa.fit_nlsLM)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 r_tref 0.150 0.00576 26.1 1.44e-40
## 2 e 0.572 0.0621 9.21 3.76e-14
## 3 eh 3.30 0.447 7.38 1.41e-10
## 4 th 36.8 0.619 59.5 2.27e-67#GRANDIS
#calculate all the TPC parameters
grandis.extra_params <- calc_params(grandis.fit_nlsLM) %>%
pivot_longer(everything(), names_to = 'param', values_to = 'estimate')
#calculate CIs for all the TPC parameters
grandis.ci_extra_params <- Boot(grandis.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(grandis.fit_nlsLM)), R = 200, method = 'case') %>%
confint(., method = 'perc') %>%
as.data.frame() %>%
rename(conf_lower = 1, conf_upper = 2) %>%
rownames_to_column(., var = 'param') %>%
mutate(method = 'case bootstrap')
#join the parameters and CIs
grandis.ci_extra_params<- left_join(grandis.ci_extra_params, grandis.extra_params)
grandis.ci_extra_params$Treatment <- "P. grandis"#meandrina
#calculate all the TPC parameters
meandrina.extra_params <- calc_params(meandrina.fit_nlsLM) %>%
pivot_longer(everything(), names_to = 'param', values_to = 'estimate')
#calculate CIs for all the TPC parameters
meandrina.ci_extra_params <- Boot(meandrina.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(meandrina.fit_nlsLM)), R = 200, method = 'case') %>%
confint(., method = 'perc') %>%
as.data.frame() %>%
rename(conf_lower = 1, conf_upper = 2) %>%
rownames_to_column(., var = 'param') %>%
mutate(method = 'case bootstrap')
#join the parameters and CIs
meandrina.ci_extra_params<- left_join(meandrina.ci_extra_params, meandrina.extra_params)
meandrina.ci_extra_params$Treatment <- "P. meandrina"#tuahiniensis
#calculate all the TPC parameters
tuahiniensis.extra_params <- calc_params(tuahiniensis.fit_nlsLM) %>%
pivot_longer(everything(), names_to = 'param', values_to = 'estimate')
#calculate CIs for all the TPC parameters
tuahiniensis.ci_extra_params <- Boot(tuahiniensis.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(tuahiniensis.fit_nlsLM)), R = 200, method = 'case') %>%
confint(., method = 'perc') %>%
as.data.frame() %>%
rename(conf_lower = 1, conf_upper = 2) %>%
rownames_to_column(., var = 'param') %>%
mutate(method = 'case bootstrap')
#join the parameters and CIs
tuahiniensis.ci_extra_params<- left_join(tuahiniensis.ci_extra_params, tuahiniensis.extra_params)
tuahiniensis.ci_extra_params$Treatment <- "P. tuahiniensis"#verrucosa
#calculate all the TPC parameters
verrucosa.extra_params <- calc_params(verrucosa.fit_nlsLM) %>%
pivot_longer(everything(), names_to = 'param', values_to = 'estimate')
#calculate CIs for all the TPC parameters
verrucosa.ci_extra_params <- Boot(verrucosa.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(verrucosa.fit_nlsLM)), R = 200, method = 'case') %>%
confint(., method = 'perc') %>%
as.data.frame() %>%
rename(conf_lower = 1, conf_upper = 2) %>%
rownames_to_column(., var = 'param') %>%
mutate(method = 'case bootstrap')
#join the parameters and CIs
verrucosa.ci_extra_params<- left_join(verrucosa.ci_extra_params, verrucosa.extra_params)
verrucosa.ci_extra_params$Treatment <- "P. verrucosa"#Join Species estimates and CIs
All_params <- rbind(grandis.ci_extra_params, meandrina.ci_extra_params, tuahiniensis.ci_extra_params, verrucosa.ci_extra_params)
All_params <- All_params %>%
mutate_if(is.numeric, round, 2)
#Plot all of the estimates
estimate.plots <- ggplot(All_params, aes(Treatment, estimate, color=Treatment)) +
geom_point(size = 2) +
scale_color_manual(name="Treatment", values=c("darkgreen","orange", "purple", "pink"))+
geom_linerange(aes(ymin = conf_lower, ymax = conf_upper)) +
theme_bw() +
facet_wrap(~param, scales = 'free_y') +
scale_x_discrete('')
estimate.plots
# #filter to only the most relavent and well characterized parameters
All_params <- All_params %>%
filter(!param=="ctmin") %>%
# filter(!param=="ctmax") %>%
# filter(!param=="eh") %>%
filter(!param=="rmax") %>%
filter(!param=="skewness") %>%
# filter(!param=="topt") %>%
filter(!param=="thermal_tolerance") #%>%
# filter(!param=="q10") %>%
# filter(!param=="e") %>%
# filter(!param=="breadth")%>%
# filter(!param=="thermal_safety_margin")
#view estimate plots
estimate.plots <- ggplot(All_params, aes(Treatment, estimate, color=Treatment)) +
geom_point(size = 2) +
scale_color_manual(name="Treatment", values=c("darkgreen","orange", "purple", "pink"))+
geom_linerange(aes(ymin = conf_lower, ymax = conf_upper)) +
theme_bw() +
labs(y = NULL)+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
legend.position = "bottom",
strip.background = element_blank(),
strip.placement = "outside") +
facet_wrap(~param, scales = 'free_y', nrow=2)+
#labeller = as_labeller(c(e = "e (Energy)", eh = " eh (Energy)", rmax= "Rmax (~nmol~O[2] ~larva^{-1}~min^{-1})",topt="Topt (Temperature °C)")), strip.position = "left") +
scale_x_discrete('')
estimate.plots
ggsave("output/TPC_estimates_SharpSchool_Respiration.pdf", estimate.plots, dpi=300, w=6, h=2, units="in")#Plot Curve and Estimate Output
#generate a combined figure of TPCs and estimate plots
figure <- ggarrange(TPC.plot , estimate.plots,
labels = c("A", "B"),
ncol = 1, nrow = 2,
heights=c(1,0.5))
figure
ggsave("output/Respiration_TPC_and_estimates.pdf", figure, dpi=300, w=6, h=8, units="in")